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Using the Liouville space framework developed in nonlinear optics we calculate the linear response 
" functions and susceptibilities of Bose-Einstein condensates (BEC) subject to an arbitrary mechanical 

force. Distinct signatures of the dynamics of finite temperature BEC are obtained by solving the 
■ Hartree-Fock-Bogoliubov theory. Numerical simulations of the position dependent linear response 

' functions of one dimensional trapped BEC in the time and the frequency domains are presented. 
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Since the first experimental observation of atomic Bose-Einstein Condensates (BEC), extensive theoretical and 
experimental effort has focused on understanding their properties [1]. BEC presents many possible avenues of research, 
as it may be viewed in a variety of ways such as mesoscopic "superatom" , gaseous "superfiuid" , or as a source of 
coherent atoms used in an Atom Laser. In particular, the close analogy between the nonlinear interaction of BEC 
matter- waves and photon waves in nonlinear optics gave rise to fascinating atom optics applications [2] . 

Optical spectroscopy has been a major tool for studying the properties of matter since the early days of quatum 
£j ' physics. And, with the advent of the laser, nonlinear optical spectroscopy has become an important technique for 
studying the properties of matter that are not accessible with incoherent light sources. The primary theoretical tool 
used in nonlinear optical spectroscopy to analyze the structure and dynamic processes in many body quantum systems 
is the linear and higher order optical response functions [3]. In this paper we extend the systematic formalism of 
optical response functions to probe the response of trapped, atomic BEC to an arbitrary mechanical force coupled 
to the atomic density. We calculate the first order (linear) suseptibilities for the condensate, non-condensate density, 
and non-condensate correlation in the time and the frequency domains. 

A major difference that arises when the formalism of nonlinear spectroscopy is applied to the mesoscopic BEC is 
that it is generally not possible to apply the dipole approximation commonly made in the atom-light interaction. It 
should also be noted that there have been a number of calculations of the linear response of BEC in the past [4]; 
the response functions calculated here provide a more fundamental Green's function description of the response to 
a force with arbitrary spatial and temporal profiles. The time domain response functions or their frequency domain 
counterparts, the susceptibilities, provide unique signatures of the dynamics of the system under consideration. 

The dynamics of zero temperature atomic BEC is commonly described using the time-dependent and time- 
independent Gross-Pitaevskii Equation (GPE). Numerical solutions of the CPE for various properties of zero tem- 
perature BEC as well as for atom optics and four wave mixing applications have been reported [5-8]. BEC at finite 
temperatures, on the other hand, require sophisticated theories that go beyond the GPE, and various approaches in- 
cluding the time-dependent Bogoliubov-de Gennes equations [13], the Hartree-Fock-Bogoliubov (HFB) theory [9-12], 
Quantum Kinetic Theory [14-19], and Stochastic methods [21-25] have been employed. 

In order to go beyond the GPE, two or more particle correlations must be taken into account. One of the original 
contributions was made by Bogoliubov with his introduction of the Bogoliubov transformation [26,27] which shows 
how the condensed state of an interacting homogeneous gas differs from that of the noninteracting gas. That result 
was extended to inhomogeneous gases by de Gennes [27]. Recently, a time-dependent version of the Bogoliubov-de 
Gennes equations was employed by Castin and Dum [13] to describe the dynamics of BEC in time dependent traps. 
Their approach was based on an expansion of the evolution equations for the atomic field operator which is valid if 
the number of noncondensate particles is small. Their result has been used for analyzing the stability and depletion 
of a strongly driven BEC [13]. 

The time- independent HFB theory (TIHFB) has been used to calculate some of the important equilibrium properties 
such as the quasiparticle excitation frequencies and the equilibrium condensate and non-condensate density profiles 
[11]. The dynamical properties of finite temperature BEC predicted by the TDHFB [12] have not been studied 
extensively, because the full solution to these equations is computationally prohibitive. The response functions 
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computed here offer a systematic perturbative approach for exploring the TDHFB dynamics, at an affordable numerical 
cost. 

The quantum kinetic theory of dilute interacting Bose gas was derived long ago by Kadanoff and Baym [14] in 
terms of nonequilibrium real-time Green's functions that parametrize the condensating gas. Their equations were 
later generalized by Hohenberg and Martin to include condensates [15] . A more contemporary version of the quantum 
kinetic theory applicable to the experimentally produced BEC has been developed in a scries of papers of Gardiner 
and Zoller [16]. They describe a system composed of interacting a condensate and a noncondensate vapor, where the 
vapor is described by a quantum kinetic master equation, equivalent to a quantum Boltzmann equation of the Uehling- 
Uhlenbeck form [20] . The master equation which describes the transfer of energy and particles between the vapor and 
condensate has been used to describe the formation of BEC. A similar formalism was also put forward by Walser et 
al. [18] who derived, using a method reminiscent of the classical Bogoliubov-Born-Green-Kirkwood-Yvon (BBYGK) 
technique, a generalized kinetic theory for the coarse-grained Markovian many-body density operator. Their result 
incorporates second order collisional processes, and describes irreversible evolution of a condensed bosonic gas of 
atoms towards thermal equilibrium. 

Important experimental observations such as the damping of elementary excitations has prompted attempts to 
develop a theory that works in the hydrodynamic regime. These have led to various versions of Quantum Kinetic 
theory such as the two fluid hydrodynamic description of finite temperature BEC [17] developed by Zaremba, Nikuni, 
and Griffin (ZNG). ZNG is a semiclassical Hartree Fock description where one neglects the off-diagonal distribution 
functions. It combines a non-Hermitian generalized Gross-Pitaevskii Equation for the condensate wave function and 
the Boltzmann kinetic equation for the noncondensate phase density. This theory was recently shown to properly 
describe damping [19]. 

It will be helpful to clarify how the various kinetic theories are related. The Green's function approach of Kadanoff 
and Baym [14] are equivalent to the equations of Walser et al. [18]. From the equations of Walser et al., the ZNG 
hydrodynamic equations may be obtained by neglecting the anomalous fluctuations. This greatly simplifies the 
equations and facilitates the numerical solution. The TDHFB equations are obtained by dropping the second order 
collisional terms from the kinetic equations of Walser et al. while keeping the anomalous fluctuations. The theory of 
Gardiner and Zoller [16] is based on the quantum Boltzmann master equations reminiscent of the quantum stochastic 
methods used in Quantum Optics; their quantum kinetic master equations also contain the higher order collisional 
processes described by Walser et al. and ZNG. A number of approximations are common to these kinetic theories: 
the Born and Markov approximation, ergodic assumption, and Gaussian initial reference distribution. The validity of 
these assumptions will ultimately be confirmed by experiments; so far the experiments have supported the predictions 
of these equations. 

The stochastic method is another approach used for the description of finite temperature BEC beyond GPE. A 
description of a Bose gas in thermal equilibrium has been developed using the quantum Monte Carlo techniques, 
based on Fcynman path integral formulation of quantum mechanics [21,22]. A stochastic scheme corresponding to 
the dynamical evolution of density operator in the positive P representation has been applied to BEC [23,24]. An 
alternative approach that describes the dynamics of the gas based on a stochastic evolution of Hartree states and 
avoids some of the instability problems of earlier works was proposed recently [25]. The stochastic methods make 
the simulation of the exact dynamics of N boson system numerically feasible; this, at present, requires assuming an 
initial state such as Hartree Fock state which is not very realistic. 

The TDHFB theory is a self-consistent theory of BEC in the collisionless regime that progresses logically from the 
Gross-Pitaevskii Equation by taking into account higher order correlations of noncondensate operators. Although 
TDHFB does not take into account higher order correlations that are done in the various quantum kinetic theories, 
the TDHFB equations are valid at temperatures near zero, even down to the zero temperature limit, and are far 
simpler than the kinetic equations which can only be solved using approximations such as ZNG. Another attractive 
feature of TDHFB from a purely pragmatic point of view is that the Fermionic version of the theory has already been 
well-developed in Nuclear Physics [9] . We therefore work at the TDHFB level in this paper and our approach draws 
upon the analogy with the time-dependent Hartree-Fock (TDHF) formalism developed for nonlinear optical response 
of many electron systems [28] . 

The paper is organized as follows: in Section II we show how to systematically solve the TDHFB equations for 
externally driven finite temperature BEC by order by order expansion of the dynamical variables in the external field. 
In Section III we define the n'th order response function in time and frequency domains and calculate the linear 
(n = 1) response function. Numerical results for the linear response functions and susceptibilities of a 2000 atom 
condensate in a one dimensional harmonic trap, and its variation with position, time and frequency, are presented in 
Section IV at zero and finite temperatures. Our main findings are finally summarized in Section V. 
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II. THE TDHFB EQUATIONS 



A. Equations of motion 



Our theory starts with a time-dependent many-body second quantized Hamiltonian describing a system of externally 
driven, trapped, structureless bosons with pairwise interactions. Introducing the boson operators a\ and a, that 
respectively create and annihilate a particle from a basis state % with wave functions 0i(r), the Hamiltonian is: 

H = H + H'(t) (1) 

where 

Ha = Y^ ( H iJ - M) 4^ + \ Yl V ^kia\a]a k ai. (2) 

ij ijkm 

The matrix elements of the single particle Hamiltonian are given by 



Hij J 



d 3 r #(r) 



0i(r), (3) 



where Vtrap(r) is the magnetic potential that confines the atoms and /j, is the chemical potential. The basis state 
</>i(r) is arbitrary; a convenient basis for trapped BEC is the eigenstates of the trap since is then diagonal. The 
symmetrized two particle interaction matrix elements in Eq. (2) are 



Vijki = - 



(ij\V\kl) + (ji\V\kl) , (4) 



where 



(ij\V\kl) = yrf 3 rd 3 r>*(r)^(r')^(r-r')^(r')^(r), (5) 
with V(r — r') being a general interatomic potential. H'(t) describes the coupling of an external field with the system: 

H'(t)= V Y i E ij (t)a\a j . (6) 

ij 

rj is a bookkeeping expansion parameter (to be set to 1 at the end of the calculation). The matrix elements Eij(t) are 
given by 

Eij = |d 3 r0*(r)^(r,t)^-(r), (7) 

where Vf (r, t) denotes a time- and position- dependent external potential that exerts an arbitrary mechanical force 
on the system. 

More general forms of H'(t) could include, for example, terms of the form y\ . Fij(t)aiaj in addition to the term 
given in Eq. (6) which couples to atomic density. In this work we shall focus on the most experimentally relevant 
perturbation; for instance, the time-dependent modulation of the trap spring constant which mimics the mechanical 
force that was recently applied experimentally [29,30]. 

The dynamics of the system is calculated by deriving equations of motion for the condensate mean field, Zi = (fii), 
the non-condensate density pij = (a\aj) — (a\)(dj), and the non-condensate correlations K%j = (aiCLj) — (a,)(aj). 
Closed equations are derived starting with the Heisenberg equations of motion for the operators, fij, a\a,j, and ajCij, 
and assuming a coherent many body state. The resulting many body hierarchy is then truncated using the generalized 
Wick's theorem for ensemble averages [31-33]: 

(Ai) * 0, (8) 
(A 1 A 2 ) = (A 1 )(A 2 ) + ((A 1 A 2 )) (9) 
(AiA 2 As) = (Ai)(A 2 )(A 3 ) + (Ai)((A 2 A 3 )) + (A 2 )((A 1 A 3 )) + (A 3 )((A 1 A 2 )), (10) 
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and similarly for products involving higher number of operators. Ai denote boson creation or annihilation operators, 
a] or Sj, and are assumed to be normally ordered. We follow the convention that normal ordered operators are time 
ordered. The double angular brackets {(AiAj)) denote irreducible two point correlations. Using this notation we have 
Pij = ((a\a,j)) and n tj = 

This procedure yields the TDHFB equations of motion [9,10,12] 

dz 

in— = [H z + f ] E(t)}z + n„z*, (ii) 

ih^ t = [h, p] - (kA* - Ak*) + r,[E(t),p], (12) 

ih lE = ^ hK + Kh *^ + ^ pA + Ap *^ + A + ( 13 ^ 

where [. . .] + denotes the anticommutator. Here, H z , H z *, h, and A are n x n matrices with n being the basis set 
size used: 

kl 

[Hz*\i,j =^2 V ijklKkl, (15) 
kl 

hij = Hij - p, + 2 ^2 Vikij \ z l z i + Pik] , (16) 

kl 

Ay = ^2 v ijki [zkzi + km] ■ (17) 
ki 

h is known as the Hartree-Fock Hamiltonian and A as the pairing field [10]. p is the chemical potential introduced in 
Eq. (2). Eqs. (12) and (13) may also be recast in a more compact matrix form [10]: 

i?i^=^[Hj,jR}+r ] E, (18) 

where H , i?, E and 7 are 2n x 2n matrices defined as follows: 

h-p A \ u _( P k \ „_/ [E(t),p] [E(t),K]+ 



H -\ a* h*-p) R -[ K * p*+i) E - {-[E(t), K ]* + -m), P ]* J- (19) 

l=(l _°, ), 7 2 = 7, (20) 



B. Solutions of the TDHFB 

Direct numerical solution of the TDHFB equations [Eqs (11 - 17)] for arbitrary field strength is complicated by 
their highly nonlinear character. A practical way to proceed is expanding all dynamical variables in powers of rj: 

z l (t) = z? ) +r l z ( P(t)+rfzf ) {i) + --. (21) 
Pi j (t)=P§ ) +VP%\t)+V 2 p% ) (*) + ■■■ (22) 

««(*) = 4 0) + + + • • • ( 23 ) 

and upon substituting these expansions in Eqs. (11 - 17), collecting all terms to a given power in 77. This results 
in a hierarchy of equations such that the equation of motion for the n'th order solution aS n \ where a denotes the 
variables z, p or k, is expressed as a function of 0, . . . (n — l)th order solutions, a^°\ , . . . , c^™ -1 ). While the zero'th 
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order equation is nonlinear, each of the higher order equations are linear; the n'th order solutions, a*-™- 1 are therefore 
obtained by solving the sequence of linear equations. This is analogous to the TDHF [28] or time dependent density 
functional algorithms used for many electron systems. 

Finding the zero'th order solution {z^°\p^°\ n^} should be the first step in solving the TDHFB. This describes the 
system at equilibrium in the absence of the external driving field, and requires the solution of the TIHFB equations: 

W<°>*<°> +W£V°>* = (24) 



- («(°)A(°)* - A«V°)*) = (25) 

(fcW/tW + «(°>/i<°>*) + (p (0) A(°) + A«V 0) *) + A<°) = (26) 

Here, , h^°\ and A' ) are nxn matrices defined in Eq. (14-17) with the variables {z,p, k} replaced by 

{z(°\p(°\ n^}. It follows from Eq. (18) that Eqs. (25) and (26) can be written in the compact form: 

7[# (0) 7,7i? (0) ] =0. (27) 

Self-consistent numerical methods for solving the TIHFB are given in the Appendix [10]. 

So far we have worked in the trap basis since this enables us to maintain a general form for the interatomic 
interactions Vijki and makes the numerical solution feasible. However, it may be more interesting to recast the 
solution in real space a^ n \r, t) (a = z, p, or k) by transforming the solution to TDHFB in the trap basis a^ n \t): 

z^\v,t) = Y,^\t)4>j{^ (28) 

3 

«W(r,t) = 2«|?)(t)^(r)^(r). (30) 

ij 

In general, real space non-condensate density and non-condensate correlations are nonlocal functions of two spatial 
points p(r',r) and n(r',r). We only computed these quantities for r = r' in this paper since these are the most 
physically accesible. Measuring these quantities with r ^ r' involves observing atomic correlations which is much 
more difficult than photon correlations. 



C. Liouville space representation 

We next introduce the Liouville space notation [3,34] that will be used in the folllowing sections. One well-known 
example for which the Liouville space formalism is used is in the optical Bloch equations, where the 2x2 density 
matrix is recast as a 4 component vector, and the Liouvillian is written accordingly as a 4 x 4 matrix superoperator. 
We rearrange the TIHFB equations Eqs. (24-26) by writing pij and Kij as vectors in Liouville space, and introduce 
the following set of n 2 x n 2 matrices (superoperators) : 

H { ~ ] -h {0) 5- -h {a) S- (31) 

rL ij,mn ~ " , im°3n n n j V 01 ) 

H {+) _ ft (0)x. +h (P)g. +V . (on) 
n -ij,mn ~ IL im u 3n ' a nj °»m ' vijmm (OZ) 

T^ij,mn — ^im^jni (33) 

^ij.mn = — A-nj ^im, (34) 

where the nxn matrices h and A were defined in Eqs. (16-17). We further define the n 2 x 1 matrices 

A ij = ^2 V mZkZi. (35) 
ki 

Using this notation, the TIHFB assume the form: 
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/ \ 




A" 




= 0. 



(36) 



Hereafter, we shall denote the zero'th order solution in Liouville space notation as ip(°\t) = 
[z<°\ z<°)* , p<°), , pW* , «(°)*] T and refer to the 2n(2n + 1) by 2n(2n + 1) matrix multiplying V^ (0) as the zero'th 
order Liouville operator, £ . 

Substituting the expansion for z i; pij, and Kij [Eqs. (21 - 23)] into Eqs. (11-13), we obtain to first order in n: 



i»^=^)(t) + C(*). 



(37) 



Here ^W(i) = p^*, « (1) *] T i.e. a 2n(2n + 1) x 1 vector with the variables in first order as its 

components. Also, 



£ = Co + d 



(38) 



where £ is the total Liouvillian, Cq was introduced in Eq. (36), and £\ and C,{t) 1 obtained by equating all first order 
terms in n are given in Appendix B. £ is given by the sum of the original TIHFB matrix £0 plus a perturbation £\. 
This perturbation induces a shift in the excitation frequencies, as will be shown below. 

Adopting Liouville space notation, the position-dependent nth order solution tp( n \r,t) can be defined using the 
relations Eq. (28 - 30) and introducing a 2n(2n + 1) x 2n(2n + 1) square matrix T(r): 



where 



fi n \r,t) = f(r)^ n >(i) 
f (r) = diag U(r), 0* (r), <D p (r), $ K (r), $;(r), $;(r) 



(39) 
(40) 



Here "diag[- • •]" denotes that T(r) is a block diagonal square matrix made of n xn blocks 4>(r), (j>*(r) and n 2 xn 2 blocks 
<I>p(r), <& K (r), $* (r), <1>* (r). <f>(r) is a diagonal matrix with the i-th diagonal element given by the basis states </^(r), and 
$ p (r) and $ K (r) are also diagonal matrices whose r/'th diagonal element are given by [^(r)]^- f ■ = </>*(r)</>j(r), and 
[$ K (r)] ij ■ ij = 4>i( r )4>i( Y ) respectively. The real space variables z^ n \r,t), p( n \r,t), and K^(r,t) are finally obtained 
by summing over the appropriate elements of the vector ip( n \r,t): 



An) 



(r,i) = E$ n) M), p(»)(r,i) 



i=i 



2n+n 2 
i=2n+l 



,(«) 



2n+2n 2 

(',*)= E ('.')■ 

i=2n+n 2 + l 



(41) 



III. THE RESPONSE FUNCTIONS 

,(n). 



The n'th order response function K K ap f(t, ti, . . . t n , r, ri, . . . r„) is defined by the relation: 

a (n) (r,t)= f [ K^(t,h,...t n ,r,r u ...r n ),V f (r 1 ,t 1 )---Vf(r n ,t n )dt 1 ---dt n dr 1 ---d 3 r n , (42) 



where a^(r, t) with a = z, p, or k are position-dependent n'th order solutions. The p subscript of Kap indicate that 
this is the response to an external field that couple to the atomic density. 

The n'th order susceptibility is defined as the Fourier transform of the response function to the frequency domain: 

poo r-oo 

X$(ft,fti,...,n n ,r,n,...r n )= / dt dh---dt n K^(t,t 1 ,...,t n ,r,r 1 ,...r n ) 

Jo Jo 

x cxp (iQt + iQih + • • • + iCt n t n ) . (43) 
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A. The time-domain linear response function 



The solution of the matrix equation Eq. (37) is 



cxp 



-j£( ( -0 



((t')dt'. 



The corresponding position-dependent solution can then be written: 



1 



dr' / T(r)W(t-t')*(r')^ (0) Vf( r/ .*')*' 



= dr' J*K^(t-t',r,r')V f (r',t')dt', 
where T(r) is as given in Eq. (40) and we have defined the 2n(2n + f ) x 2n(2n + 1) matrices 



U(t - t') = 6{t - t') exp 



$(r) = diag f$(r), $*(r), $ ( ~ } (r), $ (+) (r), $ ( - } *(r), $ (+) *(r) 



(44) 

(45) 
(46) 

(47) 
(48) 



with 6(t — t') being the Heaviside function and as discussed above, the notation "diag[- • •]" is used to denote $(r') as 
a 2n(2n + 1) x 2n(2n + 1) block diagonal square matrix with the blocks consisting ofnxn square matrices having 
the ith row and jth column given by 



[$(r)] ij= 0*(r)^-(r), 



and n 2 x n 2 square matrices 



The function K [ }\t - h, r, n) of Eq. (46), 



$(±)( r ) 



<j>*(r)(/) m (r)5 jn ± ^(r)^(r)5 im . 



^' (t - ti, r, n) = T(r)H(i - ti)$(ri)^°) , 



(49) 



(50) 



(51) 



may be viewed as the linear response function for the position dependent vector tp^ (r, t) i.e. for all the variables 
and in the trap basis. The real space response functions for the condensate z, non-condensate density 
p, and the non-condensate correlation k are therefore given by summing over appropriate indices in (t — ti,r,ri), 
using the relations Eq. (28-30): 



^)(t-ti,r,r 1 ) = £;^ ) (t-ti ) r ) r 1 ) 



i=l 
2n+n 2 



i=2n+l 
2n+2n 2 



i=2n+n 2 + l 



(52) 
(53) 
(54) 



To analyze the physical significance of the response functions it will be useful to expand them in the basis of the 
eigenvectors £„ of matrix C 



i/=l,2,...2n(2n+l). 



(55) 



We define the Green's function 
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G v {t-t') =9(t-t')exp 



rw„(t - t') 



(56) 



and further introduce r] u (r), and <L(r) as the expansion coefficients of the following vectors in the basis of 
eigenvectors £ v 



2n(2n+l) 2n(2n+l) 2n(2n+l) 

^ (0) = E *( r )£ = E ^( r )£; T(r)iL= £ <j„(r)£. 



(57) 



Wc then have 



K™(t - ii.r, n) = £ ^ r, ri)|L 



(58) 



where 



r,n) = ^ ^(r)^( ri K»G^(i-t'). 



(59) 



Here, G„(i — t') is defined in Eq. (56), ^(r), and <$„(r) arc defined as the expansion coefficients in Eq. (57), while 
the matrices $(r) and T(r) are given by Eqs. (40) and (48). Eqs. (58-59) express the linear response function Eq. 
(51) as an expansion in quasiparticle modes. 



B. The frequency domain response function 

The linear susceptibility is defined as the Fourier Transform of the response function to the frequency domain: 

poo poo 

X (1) (rj,fJi,r,ri) = / dt dt 1 K (1) (t~t 1 ,r,r 1 )eyLv{iVLt + iVL 1 t 1 ) . (60) 
Jo Jo 

It is possible to change variables t — t\ — > r and define x^ 1 )(Q,r,ri) since K^{t — ti,r,r-\) only depends on the 
time difference, t — t\. Below, we shall keep the notation K^(Sl, fii,r, ri) i.e. a function of both frequency of the 
signal, 0, and frequency of perturbation, fii, rather than K^(SI,t,ti), in line with Blocmbcrgen's notation [3]. 

In order to evaluate the Fourier Transform, we note that it follows from causality 

If 00 1 

U(t) = -— duj exp(-iut) (61) 

/GO 
dujU{oj) exp(-iujt). (62) 
-co 

Substituting Eq. (62) into the expression for K^\t — ti,r, ri) Eq. (51), and taking the Fourier Transform, we 
obtain a vector: 

^ poo roc ^ 

^(fi.fii.r.ri) = / duj dtdt 1 f(r)U(u)cxv(-iLut + iU 1 )exp{int + in 1 t 1 )$(r 1 )if {0) (63) 

* 27r * J-oo Jo 

= / (L;T(r)W(a;)$(riW 0) *(fi-wWni+w). (64) 

This implies that = and we have the susceptibility in vector form 

xf (-O i; Q u r, n) = f (r)W(ni)$(n)^ ). (65) 

The z, p and k susceptibilities are obtained by summing over the appropriate elements in the vector ( — 0; r j r i) : 
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x« si, r, ri ) = £ (-0; n, r, n) (66) 
i=i 

2n+n 2 

X W(-n;0,r,r 1 )= £ ^(-fi;n,r,r 1 ) (67) 



i=2n+l 
2n+2n 2 

X W(-n;0,r,r 1 )= £ *£(-fi;fi,r,n). (68) 

i=2n+n 2 + l 

Transforming the trap basis to the basis of eigenstates as before, one has: 

fJ ) (-fi;fl,r,r 1 )=^iC( 1 )(-fl;ar,r 1 K; (69) 

where 

4V(-n;n,r,rO = E ^' (r ^r . ( 70 ) 

and <5„(r), ^(ri), and /i„ arc as defined in Eq.(57). 

As was done in Section III A for the time domain, we have recast the linear susceptibility in two forms: matrix 
form of Eqs. (65) and an expansion in modes of Eq. (69-70). 



IV. NUMERICAL SIMULATIONS FOR CONTACT POTENTIAL 

A. Zero'th order solution (TIHFB) and the frequency shifts 

So far, all our results hold for a general pairwise interatomic interaction potential. In the following numerical 
calculations, we approximate the interatomic potential V(r — r') in Eq. (5) by a contact potential, as is standard in 
BEC applications: 

V(r-r')^U 5(r-r>), U = ^^, (71) 

m 

where a is the s-wave scattering length and m is the atomic mass. This approximation may be justified since the wave 
functions at ultracold temperatures have very long wavelengths compared to the range of interatomic potential. This 
implies that details of the interatomic potential become unimportant and the potential may be approximated by a 
contact potential. The tetradic matrices Vijki are then simply given by: 

V m = / #(r)#(r)&(r)Mr)dr. (72) 

m J J 

We assume a 2000 atom one dimensional condensate in a harmonic trap. The parameters used for our numerical 
calculations are: Uq = 47r ^ a — 0.01, and temperatures Oftwtrap/fcs and lOTk^trap/^B where w t rap is the trap frequency 
and ks is the Boltzmann constant. We used 256 grid points for position and the basis set of n — 5 states. We keep 
the trap units throughout. 

We have followed the prescription of Griffin for solving the TIHFB for ip^ in terms of the Bogoliubov-de Gennes 
Equations which are obtained from Eq. (19) by transforming to real space and using the contact interatomic potential 
[11]. In the Appendix we show the equivalence between the Bogoliubov-de Gennes Equations and the matrix H^°\ 
and summarize the numerical procedure. The eigenvalues of the non-Hermitian matrix C required for computing the 
response functions were calculated using the Arnoldi algorithm [35]. 

The eigenvalues of the Liouvillian arc trasition frequencies rather than state energies. The frequencies come in pairs 
of positive and negative frequencies; this indicates that the Liouvillian may be mapped onto a harmonic oscillator space 
for which there are always positive and negative frequency solutions. The eigenvalues of C = Co + L\ are shifted with 
respect to the corresponding eigenvalues of Co (the TIHFB equations). We list some of the representative eigenvalues, 
namely the lowest few positive eigenvalues of Co and C in Table I for both zero and nonzero temperatures. Similar 
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frequency shifts were noted by Giorgini [4]. Physically it is easy to understand how the TDHFB frequencies may be 
shifted from the TIHFB: a dynamical system contains the effect of the interacting condensate and non-condensate 
atoms which, by definition, is not present in the equilibrium system. This is analogous to optical excitations of fcrmions 
where the time-dependent Hartree-Fock theory shows excitonic shifts which are lacking by the time-independent 
Hartrcc-Fock states [28]. 

B. Time domain response 

The linear response functions were calculated by substituting the numerical solution ip^ evaluated at zero and 
finite temperatures into Eq. (51) together with Eqs. (52-54). Eq. (51) is a matrix multiplication of 2n(2n + 1) x 1 
vector ^ (0) with 2n(2n + 1) x 2n(2n+ 1) matrices T, U(t), and $ which are defined in Eqs. (48) and (40). $ and T 
are constructed in terms of the harmonic oscillator basis states which are calculated numerically from the recursive 
formula that involves the Gaussian function multiplying the Hermite polynomials [36]. The matrix U{t) was calculated 
using a MATLAB function that uses the Pade approximation for matrix exponentiation [37]. 

We first present the dependence of the linear response function on r and r' at fixed times t — t' . This gives 
a snapshot of the position-dependent correlations across the condensate. Such dependence is important since the 
experimentally produced condensates arc mesoscopic in size; in contrast, the dipole approximation usually applies in 
optical spectroscopy and consequently the spatial dependence of the response is not observable. 

In Figs 1 we display the real space response functions at times t — t' = 0/w t ra P , 7.2/cj t ra P , 15.7/cj t ra P at zero 
temperature. Fig. 2 shows the corresponding finite temperature results. Physically, t 1 and r' are respectively the time 
and postion at which the external perturbation is applied and t and r are the corresponding coordinates at which the 
measurement is made. We found that for longer times t — t' > 57r/ajtra P , the plot maintains generally the same shape 
as the third column of figures; it may therefore be possible to experimentally observe the correlations at longer times. 
At zero temperature, the correlation attains this stable shape faster than at finte temperature. 

To model a uniform perturbation applied across the condensate, we have integrated the zero temperature response 
function over r' and plotted its absolute value i.e. | J K W (r, r', t — ti)dr'\ vs. time t — t' in Fig. 3. The response 
functions K zp , K pp , K Kp grow exponentially within a relatively short time span of around 57r/w trap , so that the details 
of the structure for earlier times up to ~ 2ir / 'wtrap are not visible in the plot. This rapid growth is shown more clearly 
in the right hand column which depicts the response functions integrated over both r and r', i.e. K^\t - h). The 
figure shows the real and imaginary parts as well as the absolute values of K^(t — t{). The response function grows 
rapidly by around an order of magnitude over the time scale ~ 3tt /uJtrap- The response function keeps growing over 
time since TDHFB does not have any dissipative term; this should be a reasonable model for BEC in the collisionlcss 
regime. It also shows that, even in the absence of a dissipative term, it takes some time after the initial impulse 
at t' = before the effect of the force is reflected appreciably in the response functions. From the plots we note 
that a mechanical force applied on the condensate can be seen to "generate" noncondensate atoms and anomalous 
correlations even at zero temperature. 

The calculations of Fig. 3 are repeated at finite temperature in Fig. 4. The main difference is the smaller magnitude 
of K zp , K pp and K Kp . The fact that the BEC is less responsive at finite temperature may be attributed to the fact 
that the condensate to non-condensate interaction is greater at finite temperatures where additional collisions shield 
the effect of the applied perturbation. 

C. Frequency domain response 

The susceptibilities were computed using Eq. (65) in conjunction with Eqs. (66-68) Eq. (65) is a matrix multi- 
plication of 2n(2n + 1) x 1 vector tp^ with 2n(2n + 1) x 2n(2n + 1) matrices T, U{uj), and l>. The matrix U(uS) is 
calculated as follows: 

U{uj) = \ = Y ^ (73) 

V 

where is the right eigenvalue of L with eigenvalues io v = w„^) and £„ are the left eigenvectors such that 

J2v£v(t = 1- The eigenvalues cu u of C were calculated using the Arnoldi algorithm [35]. 

To clearly display the resonance structure, we present in Fig. 5 the absolute value of zero and finite temperature 
linear response integrated over r and ri in the frequency domain, | / Xap(—Q, r , r')drdr'\, a = z, p, k on a logarith- 
mic scale. The eigenvalues of C and hence the resonant frequencies come in positive and negative pairs; we present 
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only the positive frequencies since the function is symmetric about the zero frequency. We note that low frequency 
resonances are dominant; this may explain the absence of any oscillatory features in the time domain plots. 

Similar to the time domain calculations of Figs. 1 and 2, we show in Fig. 6 the linear response, x^( — r , r i ) , 
at zero temperature as a function of r and ri for 3 different frequencies. These frequencies represent the resonant 
frequency corresponding to the strongest peak for the condensate at O = O^trap, an off-resonant frequency at O = 
0.25aj tr ap, and the resonant frequency corresponding to the second highest peak for the condensate at f2 = 0.46oJtrap- 
In Fig. 7 we repeat the calculations for a finite temperature, and the frequencies represent the resonant frequency 
corresponding to the strongest peak for the condensate at £1 = 0.63a-> tr ap, an off-resonant frequency at fl = 1.25u>trap> 
and the resonant frequency corresponding to the second highest peak for the condensate at f2 = 1.7wtrap- Similar to 
the time domain plots, the frequency domain response is strongly position dependent. Off-resonant frequencies give 
a more complicated pattern. 

The absolute values of the zero and finite temperature X^H~ ^! r i r ') integrated over r', | j x*- 1 ^— O; f2, r, r')dr'|, 
are plotted as a function of in Figs. 8 and 9. This represents the response to a spatially uniform external 
perturbation. Since the resonance peaks vary vastly in strength we also present a contour plot in which all the 
intensities have been normalized to 1. 



V. DISCUSSION 

We have applied the systematic formalism of nonlinear spectroscopy to calculate the response functions and sus- 
ceptibilities of BEC to a mechanical force coupled to atomic density for the condensate, noncondensate density, and 
noncondensate correlations. Since our results hold for an arbitrary external perturbation, it will be interesting to 
investigate the effects of specifically taylored external force e.g. impulsive or "continuous wave" perturbations. The 
contour plots of the response functions and susceptibilities may be used for the design of experiments involving appli- 
cations of mechanical forces to a condensate. For instance, by carefully specifying the shape of the external potential 
Vf(r, t), one may effectively cancel out the response of, say, the noncondensate atoms. This would provide new insights 
into the dynamics of the condensate and noncondensate interactions. Our results show the time delay between the 
application of a mechanical force at t' = and the build up of response. We further note that the condensate and 
noncondensate atoms as well as the noncondesate correlations respond to the external mechanical force. 

The response functions computed in this paper offer a practical way to solve the TDHFB numerically with modest 
computational effort. Still, one potential bottleneck is the computational cost required to diagonalize a large matrix. 
Krylov space techniques such as the Arnoldi algorithm used in the paper considerably reduces that cost [35] . As our 
matrices were only 110 x 110 in these simulations of 1 dimensional condensate the entire set of eigenvalues could be 
calculated directly. For larger matrices, the algorithm only gives the lowest few eigenvalues. This should be sufficient 
to describe realistic experiments such as 3 dimensional asymmetric trap holding up to a million atoms. 

The damping of excitations is not included in TDHFB. In Rcf. [4] Landau and Beliaev damping were introduced 
by calculating the imaginary part of the self energy. The theory of damping of oscillations in BEC is currently far 
from conclusive and requires further investigation. This has motivated various authors to try and extend the HFB 
theory [38-40,17,19]. It should be noted that there is currently no all-encompassing theory of BEC that explains all 
observed phenomena. In addition, as noted by Leggett [41], the validity of various approximations made in some of the 
existing theories of strongly nonequilibrium dynamics of BEC (e.g. kinetics of the condensation process, the damping 
of collective excitations, and the decay of vortex states) is not entirely clear. Equations such as ZNG constitute a 
very important contribution. The TDHFB equations do constitute a systematic and consistent description of trapped 
atomic BEC at finite temperatures in the collisionless regime, just as the Gross-Pitaevskii Equation is valid near zero 
temperature and in the future it should be possible to observe directly the effects of anomalous correlations. 

Most current work on BEC excitations deals only with the linear response; however nonlinear effects should be 
observable with stronger perturbations. Just as in standard nonlinear optics, the n'th order response functions are 
expected to be most valuable for characterizing BEC in the context of matter-wave nonlinear optics. Our formalism 
allows the calculation of nonlinear susceptibilities which may be used to probe finite temperature condensates by 
four wave mixing [8]. Nonlinear response functions will be computed in a forthcoming work. Other possible future 
applications include detailed study of atom optics at finite temperatures, and superchemistry that involves the study 
of the formation of molecules from mesoscopic BEC matter- waves [42] . 
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APPENDIX A: EQUILIBRIUM TIHFB SOLUTION 



The TIHFB is given by the coupled equations: 

W<°>s< o >+W<°V (0) = 0, 



(Al) 



7[jr (0) 7,7i2 (0) ] =0, 



(A2) 



where Eq. (Al) is the time-independent GPE with and Ti^J defined respectively in Eqs. (14) and (15), and the 
matrices R (0 \ H^, 7 of Eq. (A2) are as defined in Eq. (19). 

The equilibrium solution to TIHFB is obtained variationally. Eq. (Al) is solved for z^ by using a numerical 
optimization routine to minimize the energy functional E = (H) with respect to z^* . The energy functional is found 
by applying the thermal Wick's theorem to the Hamiltonian H. For the Wick theorem, Eq. (10), to hold the state of 
the system is assumed to be described by a trial statistical density matrix of the form 



-0K 



D 



Tve-? K 



(A3) 



where the operator K is taken to be quadratic in the annihilation/creation operator for the non-condensate atoms, 
Cj/cj defined as Ci = ai — Zi/c[ = a\ — z*: 



(A4) 



with 



pf) = TvDc)c u nf> = TvDacj. 
Then one obtains for the energy E = (H): 



(A5) 



(0)* (0) . (0) 

z- z ■ + o- ■ 



ijkm 



AO)* (0)* (0) (o) 

*i *j *k rn 



Az (o)* (0) (0) 

z k Cjm 



ijkm 



, AO)* (0)* (0) (0)* (0) (0) , 2 (0) (0) (0)* (0) 
+ Z i Z ] K km "+" K i] Z k Z m + z Pik Pjra + K ij K km 



(A6) 



It is easy to show that the local minimum for E obtained by setting its derivative with respect to z^* to zero satisfies 
Eq. (Al). 

In addition, R^ is found by minimizing the thermodynamic potential for a system of bosons in thermal equilibrium. 
In Ref. [10], the generalized density matrix in thermal equilibrium, assuming a grand canonical form for the density 
matrix is shown to be given by: 



exptfHW/kT) - l 7 ' 



(A7) 



It may be shown straightforwardly that the generalized density matrix of Eq. (A7) obeys Eq. (A2), and is therefore 
a stationary solution. The proof requires using the property = 1, and the fact that a matrix A commutes with a 
function of A, f(A). The minimization of the thermodynamic potential also gives the relation [10] 



2 13 



dE 



dR. 



(o)- 



Using Eq. (A6) and the definition of R^ given in Eq. (19), Eq. (A8) implies 



hf? = -^L = Ha -n + 2j2(ik\V\lj) 



9pf 
dE 



kl 



y (0)* (0) 



.(0) 



■kl 



(A8) 

(A9) 
(A10) 
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These variational results satisfy Eqs. (16-17), derived using the Heisenberg equations of motion. 

We now discuss the self-consistent solution to TIHFB. The solution to TIHFB involves simultaneous minimization 
of Eq. (A6) coupled to Eq. (A7) where is as defined in Eqs. (16-17) and (19). Since is itself a function of 
R(°\ the solution is found iteratively. In order to evaluate R^°\ Eq. (A7), we need to diagonalizc the matrix -jH^K 
The specific form of the matrix jH^ implies that its eigenvalues and eigenvectors come in pairs. Letting V n and 
W n denote right eigenvectors belonging to the eigenvalues ±E n , 

jH^V n = E n V n , jH^W n = -E n W n , E n > 0, (All) 

where the normalisation and closure relations of the eigenvectors are: 

V n ^V m = S mn , W n ^W m = -6 mn , V n ^W m = (A12) 

and 

J2 (v n V n ^ - W n W nf ^ = 1. (A13) 

n>0 

It is noted that he eigenvectors V n and W n have the following structure: 

V n = ( yZ ) , W n =(jj2^, (A14) 

and given a right eigenvector V n , the left eigenvector belonging to the eigenvalue -E* is: V n — (U n *, — V n *), and the 
left eigenvector belonging to the eigen value —E* is: W n = (V n , —U n ). 
Using the closure relationship, of Eq. (A7) may be written as: 

R f = 12 [vrV™ + W z m W™*] n m +Y, W™W™* (A15) 

m>0 m>0 

with n m = [cxp(E m /kT) — 1]~ . The possible zero energy states are ignored, as indicated in the summation over 
m > 0. From this expression, one can calculate the matrices p and k from the explicit form of R given in Eq. (19). 
The iteration process consists of the follwoing steps [10]: 

1. Make an initial guess at Hz°^ and set p — k = 

2. Solve Eq. (Al) which determines p and z. Normalize z according to: 

N = 12 i z *i 2 +trp ( A16 ) 

i 

where N is the total number of atoms and tr denotes the trace. 

3. Solve the eigenvalue problem Eq. (All) and normalize the eigenvctors according to Eq. (A12) 

4. Calculate R from Eq. (A15) and deduce the matrices p and k. 

5. Calculate the fields TL^ and H^J and return to step 2. 

6. Stop the iteration when two successive iterations yield the same values of z, p and k to the desired accuracy. 



APPENDIX B: THE & MATRIX 



In Eq. (37), C = C + L\. The 2n(2n + 1) x 2n(2n + 1) matrix A is defined as follows: 

\ 



/ y zzL 

yzz2* 
ypzl 

~\JKZl 



yzz2 
yzzl* 
ypz2 
ynz2 



ypz2* ypzl* 
y y/tz2* yftzl* 



yzl 


y^Kh 





V z2 



o 





yzl* 



(W ph )* 

(w Kh y 



o 



(w pA y ) 



(Bl) 
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where the set ofnxn submatrices V zzl and V zz2 , n x n 2 submatrices V zl and V z2 , n 2 x n submatrices V pzl , V pz2 , 
V Kzl , and V Kz2 , and n 2 x n 2 component submatrices W ph , W pA , W Kh , and W kA of d are given as follows: 



vr = 



Vf 1 



fci — 



/ j V iklrZ k Z r V i k — ViklrZ l Z r 

kr Ir 

2j2Vukr4° ] V?%=Y,V ik lrZ r V 

r r 

V8? = 2 2 W^P*? - VrnZ^pM + E [ W?' + Virlkzf 
kr kr 

»j,fc ~~ lfeZr i ^"j v rkljZ l p ir / yrpklZj - J rVr]lkZ l 

Ir Ir 

2 ]T + V rklj zf\f + [VrikizP + V nlk z\ 



*(0) 
(0) 



y KZ } 



(0) 



.(0) 



Eh/ Jo) , t/ _(o) 



*(0) 
Prj 



Ir 



J2[v ijk izl 0) + V ijlk z\ 0) 



V KZ i 



w. 



ph 
ij,kl 



Wtki 



Ir 

2 J2 ViklJrJ - V rWP* = E KrklplT + V ^Pfr 

r r 

/ j v iklrK rj VV ij,kl — 2-^1 Vrkl 3 K ir 



E V irk,' : 



(0)* 
rj 



w. 



fcAf 
ij,kl 



VrjklK) 



(0) 



In addition, we define 



C(t) = diag 



e (-)(t) 
e(+>(t) 

V[e (+) (*)]V 



(0) 



( E{t) 

E*(t) 

e(-)(t) 

e(+)(t) 









V o o o 

We have further defined the n 2 x n 2 matrices 

£ (±) (%,fc( = S ifc (t)Jj-, ±Eij(t)S ik . 



o [*<->(*)]■ 



p(0) 



/ V« (0) */ 



(B2) 
(B3) 
(B4) 
(B5) 

(B6) 
(B7) 
(B8) 
(B9) 
(B10) 



(Bll) 



(B12) 



As mentioned in the main text, "diagLABC • • ■]" denotes block diagonal square matrix with the component matrices 
A,B,C,--- as its diagonal blocks. 



APPENDIX C: BOGOLIUBOV-DE GENNES EQUATIONS FOR CONTACT INTERATOMIC 

INTERACTION 



Griffin has provided a prescription for solving TIHFB in terms of the Bogoliubov-de Gennes Equations under the 
contact interatomic potential approximation [11]. In this section we show that self-consistent equations for [Eq. 
(A7)] is simply the Bogoliubov-de Gennes equations written in the trap basis, under the contact interatomic potential 
approximation, and summarize the numerical procedure used to find the solution to TIHFB. 

The Bogoliubov-de Gennes equations are [11]: 

H sp - p + 2U [|^(r)| 2 + n(r)] Ui (r) + U [^ 2 (r) + m(r)] w< (r) - E iUi {v) 
H sp -p + 2U [|^(r)| 2 +n(r)] Vi(r) + U ty* 2 (r) + m*(r)] ^(r) = -E iVi {v), (CI) 

where the quantities i/j g (r), n(r), and m(r) are as defined in Ref. [11] which are respectively z(r), p(r), and n(r) of 
Eqs. (28-30) when written in terms of our variables Zi, Pij, and k^. u,(r) and Vi(r) are the eigenstates to be calculated 
and can be shown to satisfy the orthogonality and symmetry relations 
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/ 



w*(r)w fc (r) - v*(r)v k (r) = 6 ik 
u*(r)v k (r) + v*(r)u k (r) = 0. 



In matrix form, Eq. (CI) is 



h A 

7 1 A * h * 



E 



(C2) 
(C3) 

(C4) 



where we have changed the basis from the position basis to the trap basis by introducing the matrix elements in terms 
of the trap eigenstates </>»(r) as follows: 



A- 



Ui = 



J 0*(r) {H» -p + 2U [|V s (r)| 2 + n(r)] } &(r)dr 
| ^(r)C/o[^(r)+m(r)]^(r)dr 

y <^*( r ) u »( r ) dr 

y 0t(r) Vi (r)dr 



(C5) 
(C6) 
(C7) 
(C8) 



Using the contact interaction, and Eqs. (28 - 30) for ip g (r), n(r), and m(r), it is clear that /i and A coincide 
with those of Eqs. (16-17); the Bogoliubov-de Gennes equations in the trap basis therefore give the same eigenvalue 
problem as that of diagonalizing the matrix 7 if of Eq. ( A7) . 

The steps to follow in solving TIHFB are therforc [11]: 

1. Solve Eq. (Al) for z\ assuming pij — = 0. 

2. Diagonalizc H of Eq. (A2) with the current value of Zi, Kij Get eigenvectors U and V. 

3. Calculate new and using U and V: 



p#0 



where N p — [exp(hLu p / kT) — 

4. Solve Eq. (Al) for z\ using the calculated values of pij and n^. 

5. Iterate: go back to Step 2. 

6. Stop the iteration when the solutions Zi, Pij and Kij converge. 



(C9) 
(CIO) 



TIHFB (T = 0hu>/k B ) 


TDHFB (T = 0hu/k B ) 


TIHFB (T = lQ%Lo/k B ) 


TDHFB (T= Whuj/kg) 























0.629 


0.3477 


0.4555 


0.3739 


0.6485 


0.3528 


0.7088 


0.6607 


0.6668 


0.7104 


0.7136 


0.6612 


0.8213 


0.7110 


1.0159 


0.8733 


0.9282 


0.8521 


1.1046 


0.8873 


1.0074 


1.0161 


1.1584 


0.8908 


1.0114 


1.0163 


1.5040 


0.9903 


1.0828 


1.1181 


1.5864 


1.0070 


1.5068 


1.1324 


1.7735 


1.0071 


1.6970 
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1.4644 


1.7760 


1.5478 


1.7666 


1.4840 


1.8211 


1.5520 


1.7788 


1.7743 


1.8832 


1.7654 


2.0183 


1.7749 


2.1995 


1.7660 


2.4251 


1.8281 


2.2553 


1.8741 


2.4716 




O /lO/IQ 

Z.4o4o 


1 QQOA 


Z.OUoo 


1.9372 


2.4935 


1.9575 


2.6454 


2.1755 


2.6050 


2.4265 


2.7705 


2.1943 


2.7907 


2.4272 


2.7927 


2.4851 


2.7914 


2.5353 


2.8689 



TABLE I. The lowest positive eigenvalues of £q and C that correspond to the TIHFB and TDHFB respectively at zero and 
finite temperatures. The eigenvalues are given in units of trap energy, TkJtrap 
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FIG. 1. K^(t — ti,r,ri) for zero temperature condensate at times t — t' = 0/^trap, 7.2/a;trap, 15.7/utrap given in different columns. 
The top, middle, and bottom rows give the response function for the condensate, non-condensate density, and non-condensate correlation 
as indicated. The dashed circle represents the spatial extent of the trapped BEC. The positions x and x' are given in harmonic oscillator 
length units. 

FIG. 2. Same as in Fig. 1, but at finite temperature 10hu>/kg. 

FIG. 3. The left column shows K^(t — ti,r) i.e. linear response in time at zero temperature integrated over ri and plotted as a 
function of t — ti and r for zero temperature. The details for the time of evolution t — t' = U7r/u;trap to 57r/u>trap are shown. The right 
column shows linear response integrated over both r and ri, (t — 1\ ), plotted as a function of t — 1\ from to 57r/w trap ; the solid, dashed 
and dotted lines represent the absolute value, real part and imaginary part respectively of the integrated response function K (1 \t - ti). 
The position x is given in harmonic oscillator length units. 



FIG. 4. Same as in Fig. 3, but at finite temperature lOTiui/ks- 

FIG. 5. Logarithm of linear response functions in frequency integrated over both r and ri vs. frequency. Top three panels - zero 
temperature, bottom three panels - finite temperature lOhcd/kg. The frequency fii is given in units of trap frequency. 

FIG. 6. X^(~ f^j r i r i) at ze ro temperature for frequencies CI = 0u>trap, 0.25u>trap, 0.46a)trap given in different columns. These 
frequencies represent the resonant frequency corresponding to the strongest peak for the condensate at f2 = Ot^trap; an off-resonant 
frequency Q = 0.25wtrap; and the resonant frequency corresponding to the second highest peak for the condensate at Q = 0.46cJtrap- The 
positions x and x' are given in harmonic oscillator length units. 

FIG. 7. X^(~ f2> Ct, r, ri) at finite temperature for frequencies f2 = 0.63wtrap, 1.25u>trapi 1.7oJtrap given in different columns. These 
frequencies represent the resonant frequency corresponding to the strongest peak for the condensate at Q = 0.63wtrap; an off-resonant 
frequency Q = 1.25wtrap; and the resonant frequency corresponding to the second highest peak for the condensate at S7 = 1.7wtrap- The 
positions x and x' arc given in harmonic oscillator length units. 

FIG. 8. The linear susceptibility at zero temperature integrated over ri and plotted as a function of f2i. The left column shows 
unsealed spectrum as a function of position; not all resonances are shown due to scaling; only the most dominant ones are represented. In 
the right hand column, the function has been normalized so that all the resonances have height of one. The position x is given in harmonic 
oscillator length units, while the frequency f2i are given in units of the trap frequency. 

FIG. 9. Same as in Fig. 8, but at finite temperature lOTiui/kB- 
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